% 8CPFSK tilted Viterbi Detection   状态:[Un_1 Vn]
id=sqrt(-1);
pi=3.14159;
h=1/8;                                       % 调制指数  
M=8;                                         % 进制数  
%bit_number=2000*3;
bit_number=50*3;
symbol_number=bit_number/(log2(M));          % 码元数目
sample_number=8;                             %每码元样点个数
Tb=1/25600;                                 % 码元持续时间
fs=1/Tb*sample_number;                       % 基带取样率

%----------------------------------------------------------------------------------------------------------------------
% data 产生0，1比特流
data=zeros(1,bit_number);
for i=1:bit_number
    temp=rand;
    if (temp<1/2)
        data(i)=0;
    elseif (temp<1) 
        data(i)=1;
    end
end;

%----------------------------------------------------------------------------------------------------------------------
% Ik={+1,-1,...+(M-1),-(M-1)} M=8  将比特流映射成M进制码元
Ik=zeros(1,symbol_number);
Uk=zeros(1,symbol_number);
Ik=[7 1 3 -3 5 -1 7 3 -3 5 -5 -1 7 1 -1 -5 5 -5 1 3 -1 -7 -5 1 -7 5 -5 -3 -3 7 -7 -5 5 -7 7 1 -1 3 -7 5 1 7 -1 1 -3 3 -7 1 3 -1];
%Ik1111=[7 7 7 7 7 7 3 5 3 -3 3 -3 7 -3 7 7 7 7 7 7 3 -3 3 -3 3 -3 7 -1 7 7 7 7 3 7 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1];
Uk=(Ik+(M-1))/2;  
%----------------------------------------------------------------------------------------------------------------------

%----------------------------------------------------------------------------------------------------------------------
% CPM g(t)=1/2LT*(1-cos(2pi*t/LT)--LRC 
L=1;
g=zeros(1,sample_number+1);
for i=0:L*sample_number
    t=i/fs;
    g(i+1)=1/(2*L*Tb)*(1-cos(2*pi*t/(L*Tb)));
end;

q=zeros(1,L*sample_number+1);
for i=1:L*sample_number+1
    for j=1:i
        q(i)=q(i)+g(j)/fs;
    end;
end;
%----------------------------------------------------------------------------------------------------------------------

%----------------------------------------------------------------------------------------------------------------------
% 计算载波相位  
% 模2*pi的绝对相位值（弧度）
theta=zeros(1,symbol_number*sample_number);                                           % 绝对相位
I_signal=zeros(1,symbol_number*sample_number); 
Q_signal=zeros(1,symbol_number*sample_number);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 利用当前码元时刻的相位增量求出相位累加
% 倾斜相位中W(与数据无关的项)
W=zeros(1,sample_number);
for i=1:sample_number
  W(i)=pi*h*(M-1)*i/sample_number-2*pi*h*(M-1)*q(i+1)+(M-1)*pi*h;
end

Uk0=0;
sigma_a=zeros(1,symbol_number);
fi=zeros(1,symbol_number*sample_number);
for n=1:symbol_number
    if n==1
       sigma_a(n)=0;  
    else
       sigma_a(n)=sigma_a(n-1)+Uk(n-L); 
    end
    sigma_a(n)=rem(sigma_a(n),8);
    theta_n=sigma_a(n)*2*h*pi;
    delta=zeros(1,sample_number);
    for m=1:sample_number
        delta(m)=Uk(n)*q(m+1); 
        delta(m)=delta(m)*4*pi*h;
        fi(sample_number*(n-1)+m)=delta(m)+theta_n+W(m);
    end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% t=nT时刻的绝对相位
for i=1:symbol_number
    fi_nT(i)=fi(sample_number*i)/pi;
end

I_signal=512*cos(fi);              % 基带I路信号
Q_signal=512*sin(fi);              % 基带Q路信号    
s=I_signal+Q_signal*id;

r_l=zeros(1,symbol_number*sample_number);
r_Q=zeros(1,symbol_number*sample_number);

for n=2:50;
    for k=1:sample_number
        s_T((n-1)*sample_number+k)=s((n-1)*sample_number+k)*conj(s((n-2)*sample_number+k));%*(cos(7/8*pi)+id*sin(7/8*pi));
    end
end

 figure(2);
subplot(2,1,1);
plot(real(s_T));
subplot(2,1,2);
plot(imag(s_T));

for i=1:sample_number            %由于发送码元前三个相同，1～8差分与9～16相同
    I_r(i)=fix(real(s_T(9))/512);
    Q_r(i)=fix(imag(s_T(9))/512);
end

for i=9:400
    I_r(i)=fix(real(s_T(i))/512);
    Q_r(i)=fix(imag(s_T(i))/512);
end

r=r_l+Q_r*id;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% viteri 检测
state_number=8;
delay=20;
len=50;   %length(s_noise)/sample_number;
transtate=[1 8 7 6 5 4 3 2;2 1 8 7 6 5 4 3;3 2 1 8 7 6 5 4;4 3 2 1 8 7 6 5;5 4 3 2 1 8 7 6;6 5 4 3 2 1 8 7;7 6 5 4 3 2 1 8;8 7 6 5 4 3 2 1];  
metric_state=-9999999*ones(state_number,delay);
survivor_state=zeros(state_number,delay);     %%其取值范围为1--64
decision_state=zeros(state_number,delay);
decision=zeros(1,len);
Ik_decision=zeros(1,symbol_number);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for n=1:len
  if n==1
    for k=1:sample_number
     r_T(k)=r((n-1)*sample_number+k);
    end 
  else
    for k=1:sample_number
     r_T(k)=r((n-1)*sample_number+k)%*conj(s_noise((n-2)*sample_number+k));
    end
  end
  %*******************************************************************************************************************
  if n==1
    %@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ 
    for i=1:state_number
        for j=i
            Vn=transtate(i,j)-1;
            Un=j-1;
            Vn_1=0;
            Un_1=0;
            refe_fi=zeros(1,sample_number);
            for kk=1:sample_number
                     refi(kk)=2*pi*h*(Un_1)+4*pi*h*(Un-Un_1)*q(kk+1);
                      ref_cos(kk)=fix(0.95*cos(refi(kk))*32768);
                      ref_sin(kk)=fix(0.95*sin(refi(kk))*32768);
                      refe_fi(kk)=I_r(8*(n-1)+kk)*(fix(0.95*cos(refi(kk))*32768))+Q_r(8*(n-1)+kk)*(fix(0.95*sin(refi(kk))*32768));
            end
            ii=1;
            metric=fix((refe_fi(ii)+refe_fi(ii+1)+refe_fi(ii+2)+refe_fi(ii+3)+refe_fi(ii+4)+refe_fi(ii+5)+refe_fi(ii+6)+refe_fi(ii+7))/3276800);
            if metric>metric_state(i,n);
                metric_state(i,n)=metric;
                survivor_state(i,n)=transtate(i,j);
                decision_state(i,n)=j-1;
            end
        end
    end
    %@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
  %*******************************************************************************************************************    
  elseif n<=delay
    %@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@  
    for i=1:state_number
        for j=1:8
            Vn=transtate(i,j)-1;
            Un=j-1;
            Vn_1=survivor_state(transtate(i,j),n-1)-1;
            Un_1=decision_state(transtate(i,j),n-1);
            refe_fi=zeros(1,sample_number);
            for kk=1:sample_number
                     refi(kk)=2*pi*h*(Un_1)+4*pi*h*(Un-Un_1)*q(kk+1);
                     ref_cos(kk)=fix(0.95*cos(refi(kk))*32768);
                      ref_sin(kk)=fix(0.95*sin(refi(kk))*32768);
                      refe_fi(kk)=I_r(8*(n-1)+kk)*(fix(0.95*cos(refi(kk))*32768))+Q_r(8*(n-1)+kk)*(fix(0.95*sin(refi(kk))*32768));
            end
                ii=1
                metric=fix((refe_fi(ii)+refe_fi(ii+1)+refe_fi(ii+2)+refe_fi(ii+3)+refe_fi(ii+4)+refe_fi(ii+5)+refe_fi(ii+6)+refe_fi(ii+7))/3276800);
                aaaa=metric+metric_state(transtate(i,j),n-1);
            if metric+metric_state(transtate(i,j),n-1)>metric_state(i,n)
                metric_state(i,n)=metric+metric_state(transtate(i,j),n-1);
                survivor_state(i,n)=transtate(i,j);
                decision_state(i,n)=j-1;
            end
        end
    end
    %@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
  %******************************************************************************************************************* 
  else 
      
    for i=1:state_number
        for t=1:delay-1
           metric_state(i,t)=metric_state(i,t+1);
           survivor_state(i,t)=survivor_state(i,t+1);
           decision_state(i,t)=decision_state(i,t+1);
        end
    end
    metric_state(:,delay)=-9999999*ones(state_number,1);
    survivor_state(:,delay)=zeros(state_number,1);
    decision_state(:,delay)=zeros(state_number,1);
    %@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
     for i=1:state_number
        for j=1:8
            Vn=transtate(i,j)-1;
            Un=j-1;
            Vn_1=survivor_state(transtate(i,j),delay-1)-1;
            Un_1=decision_state(transtate(i,j),delay-1);
            refe_fi=zeros(1,sample_number);
            for k=1:sample_number
              refe_fi(k)=2*pi*h*(Vn-Vn_1)+4*pi*h*(Un-Un_1)*q(k+1);
            end 
            ss=cos(refe_fi)+id*sin(refe_fi);
            metric=real(sum(r_T*conj(ss).'));
            if metric+metric_state(transtate(i,j),delay-1)>metric_state(i,delay)
                metric_state(i,delay)=metric+metric_state(transtate(i,j),delay-1);
                survivor_state(i,delay)=transtate(i,j);
                decision_state(i,delay)=j-1;
            end
        end
    end
    %@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
  end %else end 
  %*******************************************************************************************************************
    %--------------------------------------------------------------------
    if (n>=delay)&(n<len)
        max=1;
        for i=2:state_number
           if metric_state(i,delay)>metric_state(max,delay)
              max=i;
           end
        end
        survivor=max;
        for t=delay:-1:2
            survivor=survivor_state(survivor,t);
        end
        survivor_n=survivor;
        survivor_n_1=survivor_state(survivor_n,1);
        nnnn=survivor_n_1
%         for j=1:8
%             if survivor_n_1==transtate(survivor_n,j)
%                decision(n-delay+1)=j-1;    
%             end
%         end
        if(survivor_n-survivor_n_1)<0
            temp=survivor_n-survivor_n_1+8;
        else
            temp=survivor_n-survivor_n_1;
        end
        decision(n-delay+1)=rem(temp,8); 
        Ik_decision(n-delay+1)=2*decision(n-delay+1)-(M-1);
        if n==len-1                            % 判决(len-delay+1)时刻需保留(len-delay)时刻的幸存状态(即n=len-1时回溯到survivor_state第一时刻的状态值)
           survivor_lendelay_1=survivor_n; 
        end
    %--------------------------------------------------------------------
    elseif n==len
        max=1;
        for i=2:state_number
           if metric_state(i,delay)>metric_state(max,delay)
              max=i;
           end
        end
        survivor=1;
        for t=delay:-1:2
           survivor_n=survivor;
           survivor_n_1=survivor_state(survivor_n,t);
           for j=1:8
            if survivor_n_1==transtate(survivor_n,j)
               decision(n-delay+t)=j-1;    
            end
           end
           Ik_decision(n-delay+t)=2*decision(n-delay+t)-(M-1);
           survivor=survivor_state(survivor,t);
        end
        % 特殊处理 判断(len-delay+1)时刻
        t=1;
        survivor_lendelay=survivor;
        for j=1:8
            if survivor_lendelay_1==transtate(survivor_lendelay,j)
               decision(n-delay+t)=j-1;    
            end
        end 
        Ik_decision(n-delay+t)=2*decision(n-delay+t)-(M-1);
   end
   %--------------------------------------------------------------------
end   % for n=1:len end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

data_decision=zeros(1,bit_number);
for i=1:symbol_number
    if Ik_decision(i)==-7
      data_decision(3*i-2)=0;
      data_decision(3*i-1)=0;
      data_decision(3*i)=0;
    elseif Ik_decision(i)==-5
      data_decision(3*i-2)=0;
      data_decision(3*i-1)=0;
      data_decision(3*i)=1;  
    elseif Ik_decision(i)==-3
      data_decision(3*i-2)=0;
      data_decision(3*i-1)=1;
      data_decision(3*i)=0;  
    elseif Ik_decision(i)==-1
      data_decision(3*i-2)=0;
      data_decision(3*i-1)=1;
      data_decision(3*i)=1;  
    elseif Ik_decision(i)==1
      data_decision(3*i-2)=1;
      data_decision(3*i-1)=0;
      data_decision(3*i)=0;  
    elseif Ik_decision(i)==3
      data_decision(3*i-2)=1;
      data_decision(3*i-1)=0;
      data_decision(3*i)=1;
    elseif Ik_decision(i)==5
      data_decision(3*i-2)=1;
      data_decision(3*i-1)=1;
      data_decision(3*i)=0;  
    else
      data_decision(3*i-2)=1;
      data_decision(3*i-1)=1;
      data_decision(3*i)=1;  
    end
end
[symbolerror_num,symbolerror_rat]=symerr(Ik,Ik_decision)       % symbol error rate
[biterror_num,biterror_rat]=symerr(data,data_decision)         % bit error rate
